/***************************************************************************
 * \file   convdiff2.cpp
 * \par Description
 * This code solves the natural convection on a tall rectangular prism.
 * \f[ 
 * \frac{\partial T}{\partial t} + 
 * u_j \frac{\partial T}{\partial x_j} = \nabla^2 T
 * \f]. 
 * \par
 *  To compile and run this example type the next commands:
 * \par
   \verbatim
   % make
   % ./convdiff2 < input2 \endverbatim
 * 
 ***************************************************************************
 */

/*! 
 ***************************************************************************
 * \file   convdiff2.cpp
 * This code solves natural convection in 3D.
 ***************************************************************************
 * \author    Luis M. de la Cruz [ Sat May 23 12:06:36 CDT 2009 ]
 ***************************************************************************
 */
#include <iostream>
#include <cmath>

#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Schemes/Quick_CoDi.hpp"
#include "Equations/Momentum_XCoDi.hpp"
#include "Schemes/Quick_XCoDi.hpp"
#include "Equations/Momentum_YCoDi.hpp"
#include "Schemes/Quick_YCoDi.hpp"
#include "Equations/Momentum_ZCoDi.hpp"
#include "Schemes/Quick_ZCoDi.hpp"
#include "Equations/PressureCorrection.hpp"
#include "Schemes/Simplec.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;

// Scalar field in one dimension
typedef TunaArray<double,3>::huge ScalarField3D;

template<class T_e, class T_x, class T_y, class T_z, class T_p>
double SIMPLEC(T_e &, T_x &, T_y &, T_z &, T_p &, int, double);
void writeToFile_VTK(const ScalarField3D& data, int number, string name, double dx, double dy, double dz, string variable);

double errorT, errorX, errorY, errorZ;

int main() 
{
    double hot_wall, cold_wall;
    double length_x, length_y, length_z;
    int num_nodes_x, num_nodes_y, num_nodes_z;
    double tolerance, dt;    
    int max_iter,max_time_steps;
    double Prandtl, Rayleigh;    
    int frequency;

// Reading data from the input file
    std::ifstream input_cin ("inputConfort");
    input_cin >> length_x
	      >> length_y
	      >> length_z
	      >> num_nodes_x
	      >> num_nodes_y
	      >> num_nodes_z
	      >> hot_wall
	      >> cold_wall 
	      >> dt
	      >> max_time_steps
	      >> tolerance
	      >> max_iter
	      >> Prandtl
	      >> Rayleigh
	      >> frequency;

    double aspectY = length_y / length_x;
    double aspectZ = length_z / length_x;

    SparseMatrix<Diagonal<double, 3> > A(num_nodes_x,num_nodes_y,num_nodes_z); 
    ScalarField3D                      b(num_nodes_x,num_nodes_y,num_nodes_z); 

    StructuredMesh<Uniform< double, 3> > mesh(length_x, num_nodes_x,
					      length_y, num_nodes_y,
					      length_z, num_nodes_z);
    mesh.print();

    double dx = mesh.getDelta(X);
    double dy = mesh.getDelta(Y);
    double dz = mesh.getDelta(Z);
    int num_vols_x = mesh.getExtentVolumes(X);
    int num_vols_y = mesh.getExtentVolumes(Y);  
    int num_vols_z = mesh.getExtentVolumes(Z);  

    ScalarField3D T(mesh.getExtentVolumes());
    ScalarField3D p(mesh.getExtentVolumes());
    ScalarField3D pp(mesh.getExtentVolumes());
    ScalarField3D u(mesh.getExtentVolumes());      // u-velocity
    ScalarField3D v(mesh.getExtentVolumes());      // v-velocity
    ScalarField3D w(mesh.getExtentVolumes());      //n w-velocity
    ScalarField3D us(num_nodes_x, num_vols_y , num_vols_z); // u staggered
    ScalarField3D vs(num_vols_x , num_nodes_y, num_vols_z); // v staggered
    ScalarField3D ws(num_vols_x , num_vols_y , num_nodes_z); // w staggered

    ScalarField3D Tn(mesh.getExtentNodes());
    ScalarField3D pn(mesh.getExtentNodes());
    ScalarField3D un(mesh.getExtentNodes());
    ScalarField3D vn(mesh.getExtentNodes());
    ScalarField3D wn(mesh.getExtentNodes());
//
// Valor de los nodes en las fronteras.
//
    Range all = Range::all();
    T(T.lbound(firstDim), all, all) = hot_wall;   // Left
    T(T.ubound(firstDim), all, all) = cold_wall;    // Right

    /****/
    ScalarEquation<Quick_CoDi<double,3> > energy(T, A, b, mesh.getDeltas());
    energy.setDeltaTime(dt);
    energy.setNeumann(TOP_WALL);
    energy.setNeumann(BOTTOM_WALL);
    energy.setDirichlet(LEFT_WALL, hot_wall);
    energy.setDirichlet(RIGHT_WALL, cold_wall);
    energy.setNeumann(FRONT_WALL);
    energy.setNeumann(BACK_WALL);    
    energy.setUvelocity(us);
    energy.setVvelocity(vs);
    energy.setWvelocity(ws);
    energy.print();

    Momentum_XCoDi<Quick_XCoDi<double, 3> > mom_x(us, A, b, mesh.getDeltas());
    mom_x.setGamma(Prandtl);
    mom_x.setDeltaTime(dt);
    mom_x.setDirichlet(LEFT_WALL);
    mom_x.setDirichlet(RIGHT_WALL);
    mom_x.setDirichlet(TOP_WALL);
    mom_x.setDirichlet(BOTTOM_WALL);
    mom_x.setDirichlet(FRONT_WALL);
    mom_x.setDirichlet(BACK_WALL);
    mom_x.setUvelocity(us);
    mom_x.setVvelocity(vs);
    mom_x.setWvelocity(ws);
    mom_x.setPressure(p);
    mom_x.print();
    
    Momentum_YCoDi<Quick_YCoDi<double, 3> > mom_y(vs, A, b, mesh.getDeltas());
    mom_y.setGamma(Prandtl);
    mom_y.setDeltaTime(dt);
    mom_y.setDirichlet(LEFT_WALL);
    mom_y.setDirichlet(RIGHT_WALL);
    mom_y.setDirichlet(TOP_WALL);
    mom_y.setDirichlet(BOTTOM_WALL);
    mom_y.setDirichlet(FRONT_WALL);
    mom_y.setDirichlet(BACK_WALL);
    mom_y.setUvelocity(us);
    mom_y.setVvelocity(vs);
    mom_y.setWvelocity(ws);
    mom_y.setPressure(p);
    mom_y.setTemperature(T);
    mom_y.setRayleigh(Rayleigh);
    mom_y.print();

    Momentum_ZCoDi<Quick_ZCoDi<double, 3> > mom_z(ws, A, b, mesh.getDeltas());
    mom_z.setGamma(Prandtl);
    mom_z.setDeltaTime(dt);
    mom_z.setDirichlet(LEFT_WALL);
    mom_z.setDirichlet(RIGHT_WALL);
    mom_z.setDirichlet(TOP_WALL);
    mom_z.setDirichlet(BOTTOM_WALL);
    mom_z.setDirichlet(FRONT_WALL);
    mom_z.setDirichlet(BACK_WALL);
    mom_z.setUvelocity(us);
    mom_z.setVvelocity(vs);
    mom_z.setWvelocity(ws);
    mom_z.setPressure(p);
    mom_z.print();

    PressureCorrection<Simplec<double, 3> > press(pp, A, b, mesh.getDeltas());
    press.applyBounds(1, num_nodes_x-1, 1, num_nodes_y-1, 1, num_nodes_z-1);
    press.setUvelocity(us);
    press.setVvelocity(vs);
    press.setWvelocity(ws);
    press.set_du(mom_x.get_du());
    press.set_dv(mom_y.get_dv());
    press.set_dw(mom_z.get_dw());
    press.setPressure(p);
    press.print();    
    

    ofstream sorsumfile("./Data/sorsum"), 
      nusseltCfile("./Data/nuc"),
      nusseltHfile("./Data/nuh"),
      errorTfile("./Data/errorT"), errorXfile("./Data/errorX"), 
      errorYfile("./Data/errorY"), errorZfile("./Data/errorZ");
    
    double sorsum, nusseltC, nusseltH;
    int iteration = 0;
    
    NumUtils::interpolateToNodes(Tn, T);
    /***/
    NumUtils::interpolateToNodes(pn, p);
    NumUtils::interpolateToNodesUs(un, us);
    NumUtils::interpolateToNodesVs(vn, vs);
    NumUtils::interpolateToNodesVs(wn, ws);
    InOut::writeToFile_DX(Tn, iteration, "./Data/temp.", dx, dy, dz); 
    InOut::writeToFile_DX(pn, iteration, "./Data/pres.", dx, dy, dz); 
    InOut::writeToFile_DX(un, vn, wn, iteration, "./Data/velc.", dx, dy, dz);
    /***
    writeToFile_VTK(Tn, iteration, "./Data/temp", dx, dy, dz, "Temperature");
    /***/

    for(iteration = 1; iteration <= max_time_steps; ++iteration) {

        sorsum = SIMPLEC(energy, mom_x, mom_y, mom_z, press, 
			 max_iter, tolerance);

	if( iteration % frequency == 0 ) {
	    /***/
	    NumUtils::interpolateToNodes(Tn, T);
	    NumUtils::interpolateToNodes(pn, p);
	    NumUtils::interpolateToNodesUs(un, us);
	    NumUtils::interpolateToNodesVs(vn, vs);	    
	    NumUtils::interpolateToNodesVs(wn, ws);	    

	    InOut::writeToFile_DX(Tn, iteration, "./Data/temp.", dx, dy, dz); 
	    InOut::writeToFile_DX(pn, iteration, "./Data/pres.", dx, dy, dz); 
	    InOut::writeToFile_DX(un, vn, wn, iteration, "./Data/velc.", dx, dy, dz);
	    /***
	    writeToFile_VTK(Tn, iteration, "./Data/temp", dx, dy, dz, "Temperature");
	    /***
	    InOut::writeToFile_AM(Tn, iteration, "./Data/temp.", dx, dy, dz); 
	    /***/

	}

	nusseltH = NumUtils::calcNusseltVert(T, dx/2.0, dy, dz, 0);
	nusseltC = NumUtils::calcNusseltVert(T, dx/2.0, dy, dz, num_nodes_x-1);
	//	nusseltH = nusseltH / (aspectY * aspectZ);
	//	nusseltC = nusseltC / (aspectY * aspectZ);
	sorsumfile << iteration << "\t" << sorsum << endl;
	errorTfile << iteration << "\t" << errorT << endl;
	errorXfile << iteration << "\t" << errorX << endl;
	errorYfile << iteration << "\t" << errorY << endl;
	errorZfile << iteration << "\t" << errorZ << endl;
	nusseltHfile << iteration << "\t" << nusseltH << endl;
	nusseltCfile << iteration << "\t" << nusseltC << endl;
	
	cout << "=> Iter = " << iteration 
	     << " s = " << sorsum
	     << " e_T = " << errorT
	     << " e_X = " << errorX
	     << " e_Y = " << errorY
	     << " e_Z = " << errorZ
	     << " NuH = " << nusseltH
	     << " NuC = " << nusseltC << endl;
    }
    /****/
    return 0;
}


template<class T_e, class T_x, class T_y, class T_z,
	 class T_p>
double SIMPLEC(T_e &energy, T_x &mom_x, T_y &mom_y, T_z &mom_z, T_p &press, 
	       int max_iter, double tolerance) 
{
    double sorsum = 10.0, tol_simplec = 1e-02;
    int counter = 0;

    int pausa;

    while ( (sorsum > tol_simplec) && (counter < 20) ) {

	energy.calcCoefficients();
	Solver::TDMA3D(energy, tolerance, max_iter);
	errorT = energy.calcErrorL2();
	energy.update();

	mom_x.calcCoefficients();
	Solver::TDMA3D(mom_x, tolerance, max_iter);
	errorX = mom_x.calcErrorL2();
	mom_x.update();

	mom_y.calcCoefficients();
	Solver::TDMA3D(mom_y, tolerance, max_iter);
	errorY = mom_y.calcErrorL2();
	mom_y.update();

	mom_z.calcCoefficients();
	Solver::TDMA3D(mom_z, tolerance, max_iter);
	errorZ = mom_z.calcErrorL2();
	mom_z.update();

	press.calcCoefficients();
	Solver::TDMA3D(press, tolerance, max_iter);
	press.correction();	
	sorsum = fabs( press.calcSorsum() );

	++counter;
    }

    return sorsum;
    
}


void writeToFile_VTK(const ScalarField3D& data, int number, string filename, double dx, double dy, double dz, string variable)
{
    int Nvtk_x = data.ubound(firstDim) + 1;
    int Nvtk_y = data.ubound(secondDim) + 1;
    int Nvtk_z = data.ubound(thirdDim) + 1;
    int Ntotal = Nvtk_x * Nvtk_y * Nvtk_z;

    int w = 5;
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber + ".vtk";

    std::ofstream file (filename.c_str());        
    file << "# vtk DataFile Version 2.0 \n"
	 << "Scalar field written by TUNA \n"
	 << "ASCII \n"
	 << "DATASET RECTILINEAR_GRID \n" 
	 << "DIMENSIONS " 
	 << Nvtk_x << " " << Nvtk_y << " " << Nvtk_z;

    file << "\nX_COORDINATES " << Nvtk_x << " float \n";
    for(int i = 0; i < Nvtk_x; ++i) {
	file << i * dx << " ";
	if ( !(i % 6) ) file << "\n";
    }

    file << "\nY_COORDINATES " << Nvtk_y << " float \n"; 
    for(int i = 0; i < Nvtk_y; ++i) {
	file << i * dy << " ";
	if ( !(i % 6) ) file << "\n";
    }

    file << "\nZ_COORDINATES " << Nvtk_z << " float \n"; 
    for(int i = 0; i < Nvtk_z; ++i) {
	file << i * dz << " ";
	if ( !(i % 6) ) file << "\n";
    }
	
    file << "\nPOINT_DATA " << Ntotal
	 << "\nSCALARS " << variable << " float"
	 << "\nLOOKUP_TABLE default \n";
    for(int k = 0; k < Nvtk_z ; ++k) {
	for(int j = 0; j < Nvtk_y ; ++j) {
	    for(int i = 0; i < Nvtk_x; ++i) {
	      file << data(i,j,k) << " ";
		if ( !(i % 6) ) file << "\n";
	    }
	}
    }
        
    file.close();

}


